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Abstract. 

Continuous phase transitions occur in a wide range of physical systems, and provide a 
context for the study of non-equilibrium dynamics and the formation of topological defects. 
The Kibble-Zurek (KZ) mechanism predicts the scaling of the resulting density of defects as 
a function of the quench rate through a critical point, and this can provide an estimate of the 
critical exponents of a phase transition. In this work we extend our previous study of the 
miscible-immiscible phase transition of a binary Bose-Einstein condensate (BEC) composed 
of two hyperfine states in which the spin dynamics are confined to one dimension [J. Sabbatini 
et at, Phys. Rev. Lett. 107, 230402 (2011)]. The transition is engineered by controlling 
a Hamiltonian quench of the coupling amplitude of the two hyperfine states, and results 
in the formation of a random pattern of spatial domains. Using the numerical truncated 
Wigner phase space method, we show that in a ring BEC the number of domains formed 
in the phase transitions scales as predicted by the KZ theory. We also consider the same 
experiment performed with a harmonically trapped BEC, and investigate how the density 
inhomogeneity modifies the dynamics of the phase transition and the KZ scaling law for 
the number of domains. We then make use of the symmetry between inhomogeneous phase 
transitions in anisotropic systems, and an inhomogeneous quench in a homogeneous system, 
to engineer coupling quenches that allow us to quantify several aspects of inhomogeneous 
phase transitions. In particular, we quantify the effect of causality in the propagation of the 
phase transition front on the resulting formation of domain walls, and find indications that the 
density of defects is determined during the impulse to adiabatic transition after the crossing of 
the critical point. 



PACS numbers: 03.75.Mn, 03.75.Lm, 05.7().Fh 
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1. Introduction 

The study of non-equilibrium quantum systems is exemplified by systems undergoing a 
quantum phase transition, in which the ground state undergoes a macroscopic change of 
symmetry as a Hamiltonian parameter is changed through a critical value (TJ. These phase 
transitions often result in the formation of topological defects, occurring when the topology 
of the degenerate ground states of the new phase is non-trivial [2|. In this case, particular 
classes of excited states of the system, the topological defects, cannot be continuously 
transformed to the ground state and their existence is therefore protected by topological 
conservation laws. We can potentially unveil some of the mysteries behind out-of-equilibrium 
systems by studying the mechanism of formation of topological defects. These are imprinted 
in the system during the "turbulent" era of a phase transition but survive its subsequent 
evolution. Like fossils, topological defects provide information about a period that cannot 
be easily modelled analytically. The ability to engineer a quantum phase transition and 
recover information from topological defects is an invaluable resource for the study of non- 
equilibrium systems (3] |4] . 

The Kibble-Zurek mechanism is a theory developed by Kibble [5| and Zurek |6| with 
the goal of predicting the scaling of the density of topological defects formed with the 
rate at which the critical point of a phase transition is crossed. The theory is based on 
the assumption that every system undergoing a phase transition experiences a phase of 
impulsive behaviour where the topological defects are "crystallised" in the system. The 
breakdown of adiabaticity, due to the divergence of the response time near a critical point, 
is a characteristic of every continuous phase transition and the theory offers a model for the 
formation of defects across a wide range of systems, from cosmology [7] to condensed matter 
El|9][T0l[ni[l2l[l3l[l4l[T5][l6l. While initially applied to classical phase transitions occurring 
at finite temperature J5j|6l, more recently the scenario has been formulated and successfully 
applied to quantum phase transitions at zero temperature ifTTl [T8l [T9l l20l |2T1 l22l |23~1 l24l . 

The development of isolated, yet highly tunable ultra-cold gas experiments have provided 
exciting opportunities for the study of quantum matter and quantum dynamics. Recent 
advances in the production and manipulation of Bose-Einstein condensates (BEC) offer the 
possibility of observing a variety of quantum phase transitions, from the study of the Ising 
model with ultra-cold gases [25 1 to the exploration of the Mott insulator to superfluid phase 
transition |26|. Experimentalists in degenerate quantum gases have access to a large set of 
tools like Feshbach resonances [27], arbitrary trapping [28 1 and single-atom imaging [29 1 that 
offer them the ability to engineer Hamiltonians and to control their parameters. Furthermore 
the production of the first optically trapped BEC 11301 led to the exploration of the richer 
physics of multi-component condensates. Our ability to accurately manipulate and prepare a 
wide range of states of condensates with internal degrees of freedom allowed the observation 
of spin changing dynamics OTI and gave us access to the phase diagram of spinor BECs [ 32 1 . 
Multi-component BECs were also the key piece for the first observation of spin-orbit coupling 
in neutral atoms [33 1. 

In this work we extend a proposal for an experiment to test the KZ theory using a binary 
BEC consisting of a two hyperfine states of a single atomic species [34|. Binary condensates 
can be classified as miscible or immiscible, depending on whether the two states can naturally 
coexist in space 135). However, it has been shown that introducing a coupling between the two 
states can transform an immiscible condensate into a miscible one, and that there is a quantum 
phase transition between the two states [ 36 , 37 1 . In the strong coupling regime the mean-field 
ground state of the system consists of a superposition of the two states whose density then 
resembles the miscible phase. Here we consider a naturally immiscible condensate beginning 
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in the ground state of the strong coupling regime, which is subsequently driven back to 
the immiscible phase by quenching the coupling to zero (see Fig.[T| 1 34]. If the quench is 
non-adiabatic, regions of the system that are not causally connected will undergo the phase 
transition independently and the final configuration is a random pattern of domains of the 
two components, Fig. |TJc). The KZ mechanism is demonstrated by counting the number of 
domains formed as a function of the characteristic quench time. When applied to a BEC 
in a ring trap [28, 38, 39] this scheme is an experimentally feasible candidate to provide a 
definitive demonstration of the KZ mechanism (34). In our study we perform our simulations 
using the truncated Wigner approximation (TWA) BOl to numerically simulate the quantum 
dynamics of the BEC during the quench. This phase-space method maps quantum operators to 
stochastic variables |40| to include quantum corrections to the mean-field dynamics, allowing 
us to capture the symmetry breaking that occurs in the phase transition. 

In the vast majority of theoretical and experimental works, the KZ mechanism has been 
studied for homogeneous phase transitions |6] [17] [141 . i.e. phase transitions in which the 
critical point is traversed simultaneously for all spatial regions of the system. However, 
ultra-cold quantum gases are often confined to trapping potentials which lead to spatially 
varying phases and inhomogeneous phase transitions BT1 . In such situations additional 
considerations, such as the effects of causality, can change the quantitative details of the phase 
transition, and the KZ theory needs to be modified [42 , 43] HH. In the implementation we 
describe we can study the effects of inhomogeneities by engineering the phase transition in 
uniform systems. Spatially shaping the coupling quench, for example, allows the possibility 
of "simulating" an inhomogeneous phase transition in systems with otherwise straightforward 
dynamics such as a homogeneous ID ring BEC. The precise level of control required over 
the phase transition can be experimentally achieved with the use of modern experimental 
techniques that, for example, offer the possibility of shaping the spatial intensity profile of a 
laser using holography [45 , 46 1. Engineered Hamiltonian quenches are potentially a powerful 
tool to study the characteristics of phase transitions in a controlled fashion that can be applied 
to theoretical and experimental studies of other quantum phase transitions. 

This paper is structured as follows. In Sec. [2] we introduce the model of binary 
condensates and outline the principles behind the KZ mechanism. In Sec. [3] we model the 
miscible-immiscible phase transition in a ring BEC, and demonstrate that the number of 
defects scales with the quench time as predicted by the KZ theory l34l . We then consider 
the same phase transition in a harmonically trapped BEC in Sec. [4] We find that the counting 
of domains can be problematic, and we find a different scaling exponent with respect to the 
uniform case. In Sec. [5] we engineer a spatial quench of the coupling capable of "simulating" 
the phase transition of an harmonically trapped BEC in a ring geometry. This engineered 
phase transition allows us to verify the exponent observed in the previous section. In Sec. [6] 
we extend the method of spatially shaping the coupling to invert the process, i.e. we attempt 
to simulate a homogeneous phase transition in an harmonically trapped system. We quantify 
the effects of causality for an inhomogeneous phase transition by accurately identifying the 
subsonic and supersonic regimes of the phase transition in Sec. [7] In Sec.[8]we demonstrate 
that for this system the domain walls are seeded after the passing of the critical point, and 
finally in Sec. [9] we explore the experimental feasibility of our scheme and summarise our 
results. 
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Figure 1. Density of the two components of a binary BEC. (a) Natural ground state of a binary 
immiscible system, (b) In the strong coupling regime state (a) becomes miscible. Quantum 
noise has been added that results in the formation of domains (see text), (c) Quenching the 
coupling Q(t) to zero brings the system back to its immiscible phase. If the quench is non- 
adiabatic a random pattern of domains is created and the system is left in an excited state. We 
propose to test the Kibble-Zurek mechanism by counting the number of formed domains in 
function of the quench time. 



2. Model 



2.1. Binary Bose-Einstein condensates 

Several authors have described the physics of binary BEC mixtures |[35l l47l and coupled 
binary BECs |48 49 1. We define the Hamiltonian for our coupled ID system as 

H = H + Hi + Hc, (1) 

where Hq, Hi, and Hq are the single-particle, interaction, and coupling Hamiltonians 
respectively, defined as 



H = 



H 



H c = 



2 



t=i 



n d 2 

2m dx 2 



V(x) 



(2) 



dx 



h-i>l(x)ipi(x) 



(3) 



(4) 



The Bose field operator ?/>i (x) for component i annihilates a particle at position x, and obeys 
the commutation relations [tpi(x),ipj(x')] — S(x — x')5ij. The transversal dimensions, 
assumed to be harmonically trapped with frequency have been integrated out from the 
three-dimensional Hamiltonian, yielding the ID interaction constants = 2h 2 aij /(ma±), 
where is the scattering length and aj_ = ^h/mujj_. The coupling Hamiltonian Hq 
depends on the coupling strength and on the detuning 8 between the light-field and the 
energy difference of the states. Here we assume the coupling is on resonance and set 8 = 
for the remainder of this paper. 

The nature of the ground state of a binary system is determined by the balance between 
the interaction constants [35|. If the intraspecies interactions dominate, gu ^> <?i2, the 



Dynamics of an engineered quantum phase transition in a coupled binary BEC 



5 



energy is minimised by lowering the individual densities of both components, which then 
occupy all the available volume. In this miscible phase, the two species coexist at every point. 
However if the interspecies interaction dominates 912 ^> ga, the energy of the system is 
minimised by phase separation, which minimises the contribution of the gv2.n\ {x)n%[x) term. 
We refer to this phase as the immiscible phase. By defining 

A = ^, (5) 

9l2 

then we have 1351 

A > 1 miscible phase, 

A < 1 immiscible phase. 

However, the introduction of the coupling £1 between the states can qualitatively change this 
picture. In the strong coupling regime, when is much larger than a threshold value £! CI , the 
ground state of the system approaches the superposition state — (|1) + \2))/y/2. In this 
situation the atoms of a naturally immiscible mixture then spatially coexist as for the miscible 



phase. The calculation of the value of the critical coupling £! cr can be found in Appendix A 



2.2. Kibble-Zurek mechanism 

To formulate the KZ mechanism |5, 6| for this system, we define the dimensionless control 
parameter 

n cr - m 

e(*J = 7, > (6) 

which quantifies the distance of the system from the critical point f2 cr . For a continuous phase 
transition in the thermodynamic limit, both the equilibrium correlation length £ and relaxation 
time r of the system diverge as 

m = ktr (7) 

T(t) = R§F' (8) 

where v and z are the critical exponents of the transition and define its universality class. The 
constants £0 an d T o depend on the particular characteristics of the system, such as the atomic 
species, particle density, trapping frequencies, and so on. The correlation length represents the 
average size of the regions where the system has uniform characteristics, while the relaxation 
time characterises the time needed by the system to adjust to an external change. For a system 
dynamically approaching a critical point, there exists a moment when the relaxation time 
is equal to the characteristic time scale on which the environment is changing. The KZ 
mechanism predicts that beyond this time the system is not able to follow the quench and 
remain in quasi-equilibrium - instead it undergoes impulsive behaviour which freezes the 
configuration of defects in space. We consider a quench of the coupling parameter for our 
system of the form 

= max 0, 2fl CI (l - -^J , (9) 

where tq is the quench time. The moment when the relaxation time r[e(i)] becomes larger 
than the characteristic time of the quench \e(i)/e(i) \ determines the freezing time i through 
the condition 

r[e(t)} = e(t)/e(t). (10) 
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Solving Eq. (jT0j> with the aid of Eqs. (|7]l-(|8]l for the freezing time gives 

f _ l/(l+i/z) vz/{l+vz) nn 
T Q ■ { - Li > 

The correlation length at freezing time, £ = £,q{tq/tq) u ^ 1+,jz \ fragments the system in a 
series of regions with uniform characteristics. Since topological defects can form only on the 
boundaries between these regions |50| the number of formed defects is given by 

Ai'^aY*. (12, 



£ Co W 
2.3. Simulation method 

Our goal in this paper is to perform a computational study of the dynamics of the engineered 
quantum phase transition. In order to do so, we must go beyond the mean-field approximation 
for the system. For a system with a uniform initial density, dynamically crossing the critical 
coupling 57 CI leads to a modulational instability — such that the initial state is dynamically 
unstable. In mean-field simulations such instabilities are seeded by numerical noise. 

An approximate method for simulating beyond mean-field methods for ultra-cold gases 
is the truncated Wigner approximation ||5T1 l52l l40l . Briefly, this is a phase space method 
that simulates the evolution of the Wigner function for a system by means of stochastic 
trajectories. It is approximate in that the stochastic representation neglects some terms in 
the equation of motion involving third and higher orders derivatives of the Wigner function 
with respect to the phase space variables. However, these have a small contribution for short 
times when the number of particles in the system is much larger than the number of modes 
that are simulated Il52ll40l . 

The net result is that the equations of motion for the system are simply the appropriate 
Gross-Pitaevskii equations for the coupled condensates, but the initial state must be sampled 
from the Wigner distribution. For a weakly interacting BEC at zero temperature, this 
corresponds to populating the Bogolioubov quasiparticle modes with half a particle of 
classical noise, which represents the initial quantum fluctuations. Expectation values of 
symmetrised products of equal time quantum field operators are calculated by the averages of 
the corresponding phase space variable over a number of trajectories beginning with different 
initial conditions. 

In our simulations, the initially small "quantum noise" then seeds the modulational 
instability, and this leads to macroscopically different outcomes as might be expected to be 
realised in an experiment. Indeed, it has been plausibly argued that individual trajectories 
can be roughly interpreted as being equivalent to single experimental runs fiOl . This is the 
approach that we take for the analysis in this paper. 

The equations of motion for the two components ipi(x) for i — 1, 2 are 



ih 



&t 



h 2 d 2 

' 7) o 2 v i x ) +9u\^i{%)\ +512 1 ^3-^(2;) 1 2 
Ira ox z 



tpi{x) — HClip3-i(x). (13) 



We construct the initial states for our simulation for the homogeneous system beginning from 
the even superposition ijj+(x) — [ipi(x) + tp2(x)]/V2 with initial states tpi(x) = tp2(x) = 
y/ N/2L where N is the total number of particles and L is the length of the ring. For 
the case of an elongated BEC the initial Gross-Pitaevskii solution is found by imaginary 



time propagation of Eq. (13i in the strong coupling regime. As prescribed by the TWA, 
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the initial state is then formed by generating the complex Gaussian noises rji with variance 
Vi Vj — %/2> an d forming 

Tp(x,t = 0) = ijj + (x) + Y^[mUi(x) + r)* i v*{x)], (14) 

i 

where [ui(x), Vi(x)] are the amplitudes of the i-th Bogoliubov mode, found by solving the 
Bogoliubov-de Gennes equation for a binary BEC l37ll53l . 



3. Homogeneous Bose-Einstein condensate in a ring trap 

We begin our investigation of the KZ mechanism in the engineered phase transition in a 
binary BEC by considering a system confined in a ring trap, as was originally described in 
Ref. |34| . For our simulation we choose to simulate 87 Rb in a ring of circumference L = 96 
/itm, and take the scattering lengths to be an = 022 = 012/2 = 1.325 nm. We note that 
an and 022 are close to the measured values for 87 Rb, but that a\2 is approximately half 
the naturally occurring value, making our system strongly immiscible with A = 0.25. We 
will discuss this issue further in Sec. [9] Furthermore, we take the total number of atoms 
to be N — 5 ■ 10 4 and the transverse trapping frequency cjj_ = 2ir ■ 2 kHz which sets the 
spin healing length £ s = %/ \f2mpg s to be 1.5 times the transversal size of the system — 
enough to freeze the transversal spin degrees of freedom. The spin interaction constant is 
given by g s — (gn + g^2 ~ Igvi) /2. We simulate phase transitions with quench times tq over 
three order of magnitudes, in the range [0.1, 125] ms. The large ratio between the number 
of particles A*" and the number of simulated modes M = 1024 ensures the validity of the 
truncated Wigner approximation ll52l . 

The time evolution of the condensate density is shown in Fig. |2|a) where it is possible 
to observe an example of the final pattern of domains. The number of domains formed 
at the end of the simulation is identified as the number of zero crossings of the function 
f(x) = |V'i(a;)| 2 — IV^fx)! 2 - This method of domain counting can easily be adopted 
experimentally by performing absorption imaging of the two components after Stern-Gerlach 
separation 041 1331 . The mean number of formed domains A*d versus the quench time tq is 
plotted in Fig. |2^b). We fit the power law N$ — a-ani/ T Q to the results of the simulations 
with tq > 2 ms for which we simulated 1000 trajectories to minimise statistical uncertainty. 
We find a scaling exponent of n — 0.341 ± 0.006 while a un j = 8.63 ± 0.03. The scaling 
exponent n is in good agreement with the theoretical prediction of the KZ theory n = 1/3, 



obtained from Eq. ( 12 1 with the mean-field critical exponents v = 1/2 and z = 1 as derived 
in Appendix A and in Ref. ||49l . 

The simulation data in Fig.[2|b) shows a deviation from the expected scaling behaviour 
for fast quench times tq < 2 ms. As quenches become faster (tq decreases), the KZ 
mechanism predicts a decreasing correlation length at the freezing time, and hence a smaller 
average domain size. In contrast, in our scheme the correlation length has a lower bound 
given by the spin healing length £ s , inducing the saturation of the number of domains for 
fast quenches as seen in Fig. [2jb). The value of the saturation we extract from our data is 
compatible with the predicted maximum number given by 7V™ ax L/£ s . 

However, we also observe a mean number of domains significantly higher than the 
prediction for fast quenches [dashed box in Fig. [2jb)]. We plot the time evolution of AT d 
for these quenches in Fig.[2jd). We note that for tq < 1.5 ms the number of domains is still 
decreasing at the end of the integration time. It is in principle possible to extend the integration 
time to extract the true value of after stabilisation occurs and study phase ordering lT7TTl 
but we encountered significant numerical error and the "thermalisation" of the quantum noise 
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Figure 2. Homogeneous phase transition in a ring BEC. (a) Time evolution of one component 
showing domain formation for TV = 10 5 and tq = 47 ms. The dashed (red) line represent the 
moment when the quench reaches the critical point Q.{t) = f2 cr . The solid (green) line signal 
the end of the quench Q(t) = 0. (b) Scaling of the mean number of domains in function of 
the quench time tq in log-scale. Filled symbols represent data averaged over an ensemble of 
1000 trajectories, the open symbol over 100 trajectories. Error bars derive from the standard 
deviation. Fitting over the filled symbols yield a scaling exponent of n = 0.341 ± 0.006 
and a factor a un \ = 8.63 ± 0.03. The fitted exponent n shows good agreement between the 
numerical simulations and the KZ prediction 1/3. (c) Same data in linear scale showing the 
exponential scaling, (d) Time evolution of the mean number of domains during the quench for 
the points in the dashed box. For fast quenches the number of domains is still decreasing at 
the end of our simulation as result of phase ordering. 

when pushing the integration time beyond 600 ms fl52l . The latter is a known signature of 
the invalidity of the TWA for long evolution times, and is particularly accentuated by fast 
coupling quenches. These result in the rapid growth of the population of short wavelength 
modes, leaving the system in a highly excited state that enhances the thermalisation process. 

4. Inhomogeneous Bose-Einstein condensate in a elongated harmonic trap 

We now move to consider the same experiment performed for a binary BEC in an elongated 
harmonic trap, as is common to many experimental setups for ultra-cold gases. We take 
V(x) = mus^x 1 /2 with w x = 27r • 5 Hz, with all other parameters the same as in Sec. [3] An 
example trajectory for the trapped case is shown in Fig. [3ja). 

For the homogenous BEC it was relatively straightforward to identify the location of 
domain walls. For the inhomogeneous system, near the boundary of the condensate the 
densities |?/>i(iz;)| 2 are of a comparable magnitude to the noise that is added to the initial state, 
and the density difference function f(x) = \tpi(x)\ 2 — ^(aOI 2 can develop zeros that are not 
clearly domain walls. Experimentally an equivalent issue arises due to the presence of thermal 
and instrumental noise in the imaging system as discussed in ll56ll for quantum vortices. 

To attempt to address this issue, we introduce a new parameter 7 that limits the counting 
region for domains with a particle density larger than the threshold value p cut — 7p(0). 
Figure |3jb) shows the behaviour of the average number of formed domains in function 
of the quench time for several values of 7. The scaling exponent is again determined by 
fitting a power law to the data points for quench times where we simulated 1000 trajectories. 
Figure |3jc) plots the extracted scaling exponent n as a function of 7. Here we can see that 
n shows a dependence on the counting region with values ranging from n — 0.289 ± 0.019 
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Figure 3. (a) Time evolution of the density of one component in an harmonically trapped 
BEC. Dashed (red) line indicates the moment when Q(t) = C cr (0) (see text) while for 
the solid green line f2(i) = 0. The phase transition starts at the centre of the condensate 
where the density is highest, and propagates towards the edges, (b) Mean number of defects 
formed as a function of the quench time tq for several values of the counting threshold 
7 = Pcut/l^lCO)! 2 - Rtf is the Thomas-Fermi radius at the beginning of the quench. 
The scaling exponents are obtained by fitting a power law to the data points for which we 
simulated 1000 trajectories (filled symbols), (c) Dependence of the fitting exponent on the 
counting region expressed by 7. The dashed line indicates the exponent for an homogeneous 
quench in a ring BEC. The dashed-dotted line indicates the exponent obtained by simulating 
an inhomogeneous phase transition in a ring geometry (see text). 



to n — 0.473 ± 0.012, but is a relatively constant n ~ 0.47 for 7 > 0.3, pointing to a 
systematic difference of the scaling exponent for an harmonically trapped BEC with respect 
to the homogeneous case. For large counting regions 7 > 0.25 (where domains are counted 
in the low density region of the condensate), it can be seen in Fig.|3jb) the number of domains 
for larger quench times tq seems to converge to a value larger than predicted by a scaling 
law. This effect can be traced to the miscounting of noise as domains which results in an 
overestimate of the physical number of domains N^. 

The emergence of a different scaling exponent for 7 > 0.3 in the harmonically trapped 
case can be linked to the effects of inhomogeneity and causality, and this is the major 
topic we wish to address in this paper. Qualitatively, the value of the critical coupling is 
connected to the particle density, therefore the phase transition in an harmonically trapped 
BEC becomes inhomogeneous. The centre of the BEC, the point of highest density, enters 
the new phase first and the transition proceeds then towards the edges, as can be seen in 
Fig.[3ja) where the domain formation near the edges occurs at later times, and a propagating 
phase front can clearly be identified. In a trapped system the spatial dependence of f2 cr (x) 
and consequently of the control parameter e tI (x, t) breaks the translational symmetry giving 
a preferred direction of motion for the newly formed domains. In systems of finite size this 
has the effect of increasing the annihilation rate of domains or the rate at which they escape 
our counting regions. This effect becomes increasingly important for slower quenches, and 
increases the observed scaling exponent. 

4.1. Causality in the harmonically trapped system 

The issue of causality is related to the speed of the front separating regions in the unbroken 
(old) and broken (new) symmetry phase. It has been proposed ll42l l57l 1581 that in the case 
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where the front is slower than the speed of sound, information can propagate from regions 
in the new phase to regions in the old phase. This exchange of information influences the 
choice of symmetry of the system in the old phase, reducing the number of defects that are 
eventually formed. To analyze this phenomenon, we start by deriving an expression for the 
density of a binary condensate in the strong coupling regime. Applying the Thomas-Fermi 
(TF) approximation for the density ll59l we have 

fi — V(x) + M2(0) 



Pi{x) = p 2 {x) = p(x) = 



.912 + .9 



(15) 



for x smaller than the Thomas-Fermi radius _Rtf defined by the solution of fj, — V(Rtf) + 
hfl(Q) = 0, where /i is the chemical potential resulting from N = 2 J dxp(x). As described 
above, the critical coupling £l cl is now spatially inhomogeneous, with dependence given by 
hfl cr (x) = (gi2 — g)p(x). We can modify our definition of the control parameter Eq. (j6j) to 
reflect the spatial inhomogeneity of the phase transition 



ttr(x,t) 



n CI (x) - n(t) 



1 



,P(0) 



1 



t 



(16) 



Q cr (x) p(x) 

where we have used f2(t) — max[0, 2O cr (0)(l — t/ro)]. In inhomogeneous phase transitions, 
different regions of the system experience different effective quench time as a consequence 
of the spatially dependent critical coupling. The effective quench time is again defined as the 
change rate of the newly defined control parameter e tI (x, t), 



tq(x) 



de tI {x,t) 



dt 



= T Q~ 



n - v(x) + nn(o) 



(17) 

i?TF ■ Therefore, outer 



2[fx + nn{p)} 

We note that the effective quench time tq(x) approaches zero for x 
regions of the BEC are subject to a smaller correlation length at freezing time, and the average 
size of the domains decrease in the proximity of the condensate edges, see Figs.|3ja) and[4] 

Solving the equality e t i (^f, tp) = for xp defines the trajectory of the front which in 
the case of harmonic potential moves according to 



x F (t) = 



-\P- 



nn(o)} ( — -l 



(18) 



We have verified that Eq. ( 18 1 accurately reproduces the trajectory of the front by comparing 
it with the simulations, as shown in Fig. [4] The speed of the front during the transition follows 
from Eq. ( 1 8 1 as 

2[M- 



vf(xf) 



dt 



mu^TQXp 



(19) 
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4.2. Kibble-Zurek scaling exponent for the harmonically trapped system 

In order for the new phase to transfer information to the old phase, the front speed has 
to be smaller than the relevant speed of sound at freezing time given by i = £/f = 

&/[ro (1+ " ) rQ { *~ 1) ]^3 g2| where we have used £ = £ e(*)~ v = ?o[tq/t ]tt^ and 

f = TQe(t)~ yz = [toTq 2 ] 1 +" 2 = . For the phase transition under consideration here (z = 1), 
v is independent of the quenching time and equal to the local speed of sound v s (x) — 
\/(g + gi2)p(x)/{2m) |6U]|37]|. Solving the inequality 

2[p + nn(0)} < j (g + g 12 )p(x) 

for x defines the region X within which the formation of defects is suppressed. For the 
parameters used in this work Eq. d20| has no solution for tq < Tq 150 ms. Thus for a 
large range of quench times we simulate, the front is always faster than the sound, and the 
phase transition is not affected by the broken symmetry choices in the neighbouring regions. 
For tq > Tq on the other hand, the speed of the front equals the speed of sound in two 
points (see Fig. [5^, where the pink solid line compares the speed of the phase transition front 
to the local speed of sounds for tq = 223 ms.) When the front first enters the condensate at 
the centre it moves faster than the sound for a short distance, until it slows into a subsonic 
regime. Approaching the edges, where the speed of sound goes to zero, the front again 
becomes supersonic. As the size of the region with suppressed domain formation increases 
with the increasing quench time, the effect of causality on inhomogeneous phase transitions 
is to increase the scaling exponent, by further reducing the total number of domains formed 
compared to the fully supersonic case. We can compute the number of formed domains in this 
scenario by integrating the effective quench time tq (x) over the supersonic region: 

iV d oc / dxT Q (x)- lJ /<- 1+ ^ + / dXT Q {x)- v / { - 1+vz \ (21) 

JO Jx 2 

where x\ and £2 are respectively the inner and outer boundaries of the subsonic region X. 



When integrated numerically with v = 1/2 and 2 = 1, Eq. (21 1 yields a scaling exponent 
n«0.37. 

In more detail, we note that the comparison of the local speed of sound to the front 
propagation speed does not provide a rigorous criterion for delimiting the suppression of 
domain formation. In fact, the authors of Ref. [57 , 58 ] note that defect formation can occur for 
vf/v s > 0.5. However, the introduction of any proportionality constant in Eq. |20| modifies 
the value of the critical quench time Tq, whereas the value of the scaling exponent remains 
n s» 0.37. We address the issue of the critical value for the speed of the front in Sec. [7] 

4.3. Spatial density of defects 

We turn now to the analysis of the dependence of the spatial density of defects on the 
condensate density for a harmonically trapped BEC. A comparison between the front speed 
for a range of quench times and the local speed of sound is shown in Fig.|5ja). In Fig.|5jb) we 
plot the mean spatial distribution of the defect density at the end of the integration time for the 
same quench times as in Fig.|5ja). We can see that for fast quenches (for which the front speed 
is always larger than the speed of sound) the domain density increases towards the edge of 



the condensate as suggested by Eq. ( 17 1. The rise in defect density for x > 0.87?tf is due to 
the difficulty of distinguishing noise from domains in the low density region, as described in 
Sec. [4] Figure|5jb) shows that for slower quenches, the density of domains begins to decrease 
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Figure 5. (a) Comparison between the speed of sound (dot-dashed line) and the speed of 
front propagation for the quench times expressed in (b). (b) Final spatial density of defects 
for the trapped system for a range of quench times tq as indicated in the figure. i?TF is the 
Thomas-Fermi radius of the condensate. 



for increasing x, rather than increasing as is observed for faster quenches. This is likely due 
to the suppression of domain formation as described above. The onset of the suppression 
is expected to result in a sudden and sharp decrease in the number of domains. However, 
we have been unable to observe this clearly in the simulation results, as the positions of the 
domains are not fixed when > 0, and so domains initially seeded in the supersonic region 
are able to move into the subsonic region before they become clearly distinguishable. 



5. Inhomogeneous phase transition simulated in the ring geomety 

While the simulations of the experiment in the harmonic trap suggest a scaling exponent for 
the number of defects of n ~ 1/2, this is still open to question. To confirm the modified 
scaling exponent derived in Sec. [4] we now "simulate" the inhomogeneous phase transition 
of the harmonically trapped BEC. We acheive this by performing numerical simulations of 
a binary BEC in a ring trap subject to a spatially dependent quench of the coupling strength 
designed to reproduce the physics of the trapped system. We use 



fi(x,t) 



P( x ) V T Q 



(22) 



where p(x) is the Thomas-Fermi density of the trapped BEC system we are "simulating" 
[Eq. (fl5)], and the circumference of the ring L — 110 fim < 2i?TF (equivalent to 7 = 0.15) so 
that we avoid any divergence of Cl(x, t) at the edge of the simulated system where p(x) — > 0. 
In this situation the propagating phase transition front is due to the spatially dependent 
coupling parameter Q(x,t), rather than the spatially dependent density. By design, this 



simulation has the same control parameter for the quench as et r (x, t) in Eq. ( 16 1. However, 
we avoid the domain-counting problem described in Sec. [4] as the BEC has constant density 
around the ring. 

The evolution of the density of one component of the BEC is shown in Fig. [6ja). The 
critical coupling is first reached at x — 0, and the front of the phase transition moves around 
the ring in exactly the same manner as in the trapped BEC with a spatially constant coupling. 
The scaling of defect number for this simulation is shown in Fig. [6jb), and we find a scaling 
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Figure 6. (a) Evolution of the density of one component during an inhomogeneous coupling 
quench in a ring trap with tq = 75 ms. The red dashed line shows the approximate motion 
of the front of the phase transition from the solution of ei r (x,t) = 0. (b) Scaling of the mean 
number of domains iVd as a function of the quench time tq on a log-scale. The fitting with 
a power law results in an exponent n = 0.497 ± 0.015. The inset shows the same data on a 
linear scale. 



exponent of n = 0.497 ± 0.015. Thus we conclude that the exponent of n ~ 0.5 that was 
determined in the harmonic trap is robust. 

5.7. Kibble-Zurek scaling exponent for the inhomogeneous phase transition 

We now estimate the Kibble-Zurek scaling exponent expected for the inhomogeneous phase 



transition in a ring BEC. For a ring BEC the condition Eq. (20 1 for domain suppression due 
to causality becomes 



2[p + nn(0)} (g + g 12 )N 

n < V 7, T ' (Zi > 

where on the right hand side we have inserted the constant speed of sound in the ring l37l . 
The number of domains scales as function of the quench time as 

iV d oc f dxT Q (x)-»/( 1+ ^\ (24) 



where x is the solution to the equality of Eq. (23 i. If the size of the system L is such that the 



inequality Eq. ( 23 1 is never satisfied, the temporal dependence in Eq. ( 24 1 factors out. This 
will increase the number of domains compared to the case of the homogeneous transition, 
but leaves the scaling exponent unchanged at n = 1/3. However, if there is a region where 
the speed of the phase transition front is less than the speed of sound, the defect-suppression 
region grows with the quench time and the scaling exponent is n = 4/3. 

For the results we present in Fig. |6|b) the quench times are small enough that there is 
no region of domain suppression. However, we find the scaling exponent to be n ~ 0.5, in 
disagreement with the prediction above. We believe this discrepancy is caused by the breaking 
of translational invariance, as we describe below in Sec. |5.2| 
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Figure 7. Effect of inhomogeneity on the decay of domains. (Left panels a-d) Evolution of 
the condensate density for one component in a ring trap with 16 seeded domains. Left column 
(a,c) is for the the homogeneous quench given by Eq. \26\ , and the right column (b,d) is for 
the inhomogeneous quench given by Eq. (1271. The top row (a,b) is for a fast quench tq = 15 
ms, where it can be seen that there is no change to the domain pattern that is seeded. The 
bottom row (c,d) is for a slow quench tq = 286 ms. It can be seen in (d) that there is 
significant dynamics and domain annihilation for the inhomogeneous quench, (e) The number 
of domains surviving an inhomogeneous (blue curve) and a homogeneous (red curve) quench 
in function of the quench time tq . 




5.2. Translational invariance and domain annihilation 



For a BEC in a ring trap the system is symmetric under rotations about the center of the 
ring. This is reflected in the fact that the Hamiltonian Eq. ([T]i of the main text is invariant 
under translations x — > x' + a when V(x) — (the definition of a homogeneous system). 
However, harmonically trapped condensates and spatially inhomogeneous quenches of the 
coupling break this symmetry. To clearly demonstrate the effect of breaking translational 
invariance on the phase transition we simulate both a homogeneous and an inhomogeneous 
quench in a uniform system, but seed the domain formation "by hand." 
We begin with the initial state 



[l/>l(x),-0 2 (x) 



1 + A sin A; 



2irx 



1 - A sin 



, 2irx\ 



where k is an integer that determines the number of seeded domains, and A 
amplitude of the seed. The coupling is 

t 



max 



0, u (x,t) 
for the homogeneous quench and 



0.3- 



111JLX 



0,2fi r 



p{x) 



0.3- 



t 



(25) 
0.1 is the 

(26) 
(27) 



for the inhomogeneous quench emulating the phase transition in the harmonically trapped 
BEC. Both choices ensure that £l(x, 0) < f2 cr everywhere at t = 0. The domains will grow 
from the seed with no propagating front for the phase-transition, but the nonzero coupling will 
affect the domain dynamics. 

We plot typical results in Fig. [7] For fast quenches we find that the domains do not have 
time to move or decay, and the number of domains observed at the end of the integration 
time is the same as the number seeded for both homogeneous and inhomogeneous quenches. 
However, for slow quenches we find that £lj(x, t) results in a systematically smaller number 
of domains than £ln(x,t) for the same quench time, due to the motion and subsequent 
annihilation of domain walls. This effect will increase the scaling exponent for the number of 
domains, and for our parameters it increases it to n « 0.5. 
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6. Simulating a Homogeneous Phase Transition in an Harmonically Trapped BEC 



The success of the quantum "simulation" in Sec.[5]suggests that we can potentially engineer a 
spatially-dependent coupling Q(x, t) such that the phase transition in an harmonically trapped 
BEC happens everywhere simultaneously — that is, the phase transition is homogeneous. We 
describe our efforts to realise this below. 

A harmonically trapped BEC reacts to a coupling quench by changing its shape. This 
means that a coupling quench of the form Q(x,t) = m,ax[0, 2Q® T (x)(l — t/ra)], where 
0° r (a;) is the initial critical coupling, will not produce an homogeneous quench as it fails to 



take into account the change of shape of the system density. We extend Eq. ( 15 1 to take into 
account a spatially and time dependent coupling with the following ansatz for the density 

p(t) - V(x) + hft(x,t) 



p(x,t) = 



(28) 



5i2 + .9 

This reflects Eq. ( p~5] > in the case where the coupling is spatially inhomogeneous and uses the 
local density approximation, i.e. the condensate density is influenced only by the local value 
of the trapping potential and coupling between the components. We also assume the coupling 
£l(x, t) to be proportional to the critical coupling at any time 

m(x, t) = c(t)hn cr (x) = c(t)[g 12 - g] P (x, t), (29) 

where we have used the dependence of fi cr (a?) on the density. We have also introduced the 
function c(t) = 2(1 — t/rq) which describes the proportionality between fl(x, t) and fl cr (x). 
Eqs. ( |28|j29l ) constitute a system of equations describing the mutual relation between the 
condensate density p(x, t) and the coupling Q(x,t) that we set to a time-varying multiple of 
the critical coupling fi cr (a:, t) decided by p(x, t). The condition expressed by Eq. (29 1 ensures 
the simultaneous crossing of the critical point for c(t) = 1. 

Substituting Eq. (29 1 into Eq. (28 i and solving the equation J dxp(x) = N/2 we obtain 
an equation for tt(x, ijfor an harmonically trapped system 

\ 2 /3 



nn(x,t) 



c(t)(g 12 - g) 
G{t) 



3N 



-G{t) 



2 2 

mcuzx 



(30) 



where G(t) = g\ 2 + g — c(t){g\2 — g). We have implemented the coupling quench of 
Eq. ( pO] ), and the evolution of the density of one of the components is shown in Fig. [8] We 
find that the phase transition is effectively homogeneous in the range \x\ < 50 pm. However, 



in deriving Eq. ( 30 1 we have assumed that the condensate adiabatically follows the spatial 
dynamics of the quench at any time, i.e. the quench time is much larger than the condensate 
reaction time. For faster quenches the condensate fails to follow the changes in the coupling 



and the conditions leading to Eq. ( 30 » break down. As tq decreases, the region where the 
transition is effectively homogeneous rapidly shrinks or even disappears, complicating the 
process of domain counting. This leaves us with a narrow range of quench times tq with a 
homogeneous transition, which is insufficient to reliably extract a scaling exponent. 

From Fig. [8] it is also evident that the quench of Eq. (30 1 results in the onset of breathing 
dynamics that further complicates the counting process. In order to explain the origin of the 
breathing mode we derive the Thomas-Fermi radius of the condensate in the approximation 
given by Eqs. ( |28p9| ). Inserting the solution of the system of equations into the normalisation 
condition N = 2 JJ^ F p(x)dx, we obtain the chemical potential 



(l - c(t) + \/A - c(£)a/A 



2/3 



1/3 



jV 2 / 3 



(31) 
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Figure 8. Simulation of a homogeneous phase transition in a trapped BEC. Time evolution of 
the particle density of one component for a quench time of tq 79 ms and N = 20000. The 
dashed green line denotes the region where the phase transition is homogeneous. The inset 
shows the time evolution of the coupling Q(x,t). 



and the Thomas-Fermi radius 
"3 



i?TF — 



312 (l - c(t) + a/A - c(t)y/A 



1/3 / 2 \ -1/3 



iV 1/3 . (32) 



From (32 1 it can be seen that i?TF increases as c(i) — > 1 + , explaining the breathing in the 



case of a coupling quench from c(0) > 1 to c(tfi na i) = 0. 
7. Threshold Value for Causality 

Previous analysis of the effect of causality on the formation of defects found that defect 
suppression occurs when the speed of the front is a fraction of the speed of sound If57ll58"l . 
In this section we determine the value of the critical velocity of the phase front below which 
the formation of domains is suppressed for the coupled binary BEC. To do so we simulate 
the inhomogeneous phase transition in a ring with a front moving with the velocity vp. We 
control the velocity of the front by using the following quench of the coupling 

KQ(x, t) = max 0, fl cl - 12 + tanh (x — — — vptj + tanh ^— x — — vptj j , (33) 

where the initial distance between the fronts w is introduced to avoid the effect of the front 
entering the system. We simulate this system using the same parameters as Sec. [3] with 
w = L/20 for a time of t = 150 ms. 

The mean number of defects formed as a function of the speed of the front vp is shown 
in Fig. [9] which indicates a value for the critical speed of vp/v B f=s 0.37. In the mean field 
approximation the transition between the supersonic and subsonic regime is expected to occur 
at a precise fraction of the speed of sound. However the presence of fluctuations in the density 
of the system means that there is some uncertainty in the the local speed of sound, smearing 
out the exact value of the transition. In Ref. [58] the authors determine the relevant speed for 
the onset of causality using a moving step function for the control parameter. In this work 
we use a spatially continuous control parameter and consequently the phase transition occurs 
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Figure 9. Determination of the critical speed for domain suppression. The mean number of 
domains formed is plotted as a function of the speed of the phase transition front up , resulting 
in a critical speed of vp/v B ?s 0.37. The inset shows the coupling quench used to test the 
critical velocity. 



on a finite width. This finite width is however small and comparable with the spin healing £ s 
hence it is unlikely to have substantial effects. 

8. When is the Defect Density Determined? 

Historically the KZ mechanism has been at the centre of a debate about the timing of the birth 
of defects 16T1 l62ll . On one hand, it has been suggested that the appearance of defects can 
be traced to the fluctuations that freeze out at a time t before the transition occurs, i.e. when 
the dynamics first changes from being adiabatic to impulsive [63 1. However, there is also a 
symmetric point after the transition, when the system dynamics again becomes adiabatic. An 
alternative viewpoint has been expressed that it is this point, after the transition has occurred, 
when the eventual defect density is determined lf62l l70ll . 

In our numerical simulations we can probe the time at which the domain density is 
determined by employing a piecewise linear coupling quench for the ring-shaped system, 
with a quench time tq that is different on each side of the critical point. We choose a quench 

v ; \ 2fi cr (i - */tq) if n(t) < n cr . 



For a given quench time tq, from Eq. ( [34} we can see that the rate of change of £1 before the 
critical point is 3/2 times smaller compared to the critical point. The choice of the relatively 
small factor of 3/2 is due to the finite range of quenching times we can reliably simulate. 
Very fast quenches result in a larger number of domains that relaxes even after the end of the 
quench; on the other hand, very slow quenches are more likely to experience problems with 
the truncated Wigner approximation (Sec. [3J. 

We performed simulations of this system with the parameters otherwise as described 
in Sec. [3] and in Fig. 10 we plot the number of domains formed as function of tq for this 



asymmetric quench. We fit the data with a power law iVa = a asy /rq and derive the factor 
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Figure 10. Mean number of domains formed for the asymmetric quench of Eq. (open 
red diamonds) compared to that of the symmetric quench (open blue circles, black dots) as 
in Fig. |5Jb). The two sets of data share the same quench time after the critical point. The 
agreement of the results suggests that the domains are imprinted after the critical point. The 
inset illustrates the difference in the coupling as a function of time. Solid blue line: symmetric 
quench. Dashed red line: asymmetric quench. 



Oasy = 8.74 ± 0.06 which is in good agreement with the same quantity extracted from the 
data of Sec.|3]in Fig.^c) of a uni = 8.63 ± 0.03. 

The KZ mechanism does not provide a conclusive prediction for the total number of 
topological defects formed in a phase transition, instead focusing on their scaling. This 
uncertainty is usually treated by saying that the typical distance between defects is /£, where 
/ is unknown and typically varies from 1-10 l64l . However, in simulating the quenches in 
Fig. [TO] we have only altered the slope of the ramp before the critical point, while all the other 
parameters remained identical. As a result it seems reasonable to assume that the pre-factor / 
is the same for both scenarios. 

With the assumption that / is fixed, if the defects were decided before the critical point, 
the absolute number of domains for the scenario with quench time 3/2 times slower should 
have been reduced by a factor of (3/2)™ = 1.14, where we have used n = 1/3. Then 
we would have expected to measure a prefactor of a asy w 7.54. This is well outside the 
uncertainty range of the value that we have measured. The equivalence of the total number 
of domains for both type of quenches strongly suggests the defect density is decided after the 
transition. 

This conclusion is compatible with the physical mechanism of domain formation in our 
system. In the miscible-immiscible quantum phase transition, domains form as a consequence 
of modulational instability of the system l65l . In this situation the population of unstable 
modes with complex eigenvalues grows exponentially after the transition, suggesting that the 
number of formed domains is decided by the characteristic length of the most unstable mode 
at the moment of the return of adiabaticity. 
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9. Experimental Feasibility 

We now examine the experimental feasibility of our scheme. The miscible and immiscible 
phases were experimentally observed in various binary systems |66| and recently both phases 
were realised with the same pair of atomic species using Feshbach resonances [67|. The 
results we present in this paper are obtained in the regime where the two components are 
strongly immiscible with A = 0.25. This means that the system spin healing length is 
relatively short, and leads to both a large number of domains and their straightforward 
identification. While such a strongly immiscible system may be difficult to realise in practice, 
it does not present an overwhelming obstacle to the experimental realisation of the scheme 
described here. 

No pair of hyperfine states of 87 Rb and 23 Na naturally satisfy the criteria of naturally 
strong immiscibility, but the combination of the \F = 1, m-p = +1) and \F — 2, m-p = — 1) 
hyperfine states in 87 Rb has an interspecies Feshbach resonance |68| that can be used to tune 
A to 0.8 while keeping <?n w 322- Although the use of a Feshbach resonance enhances losses, 
smaller values of A accelerate the process of domain formation allowing for the observation 
of the density pattern. We estimate that for a ring BEC with L = 50 pm and N = 5000, 
A = 0.8 yields £ s = 0.4 p,m, t = 351 ms and iVf 1 ax « 50. 

Recently, Lin et al. reported the observation of the miscible-immiscible quantum phase 
transition in spin-orbit coupled Bose-Einstein condensates l33l . The group coupled two 
Zeeman sub-levels of the F = 1 hyperfine state of 87 Rb and was able to measure the phase 
separation of the dressed states across the critical point by changing the coupling strength 
£1 The method is not affected by increased inelastic losses, as outlined in the previous 
paragraph for Feshbach resonances, and is capable of achieving higher separation regimes 
(A <C 1) that make the task of counting domains easier. Similar results have also been 
reported by Nicklas et al. (69'). However, we note that the spatial configuration of the dressed 
state is not directly accessible, but is instead reconstructed from absorption imaging of the 
bare components. Hence, the higher separation and stability come at the price of a more 
complicated detection process for the number of domains. In addition, the results presented 
in Ref. ||69l and our own numerical simulations suggest that the appearance of domains in the 
dressed state is slower than in the procedure proposed here. As a consequence, this alternative 
scheme could potentially result in higher domain loss before the counting, further altering the 
scaling exponent. 

10. Conclusions 

In conclusion, we have shown that the number of defects formed in the coupling induced 
miscible-immiscible phase transition in a ring-shaped binary BEC scales as predicted by 
the Kibble-Zurek theory. The scaling exponent we determine numerically for the number of 
defects for this system agrees with that predicted by using the mean-field critical exponents. 
The results suggest the scheme is a good candidate for the experimental testing of the 
predictions of the KZ mechanism in an experiment with ultra-cold gases, taking advantage 
of the system's isolation from the environment and the precise control of the parameter fl that 
can be provided by microwave or laser coupling. 

We have also examined the phase transition in a binary BEC in an elongated harmonic 
trap, and found that this yields a modified scaling exponent compared to the homogeneous 
case. We have shown how an engineered coupling quench can be used to simulate 
inhomogeneous phase transitions in a ring BEC, and derived the scaling exponent of the 
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harmonically trapped case independent of the domain counting problem inherent to trapped 
BEC systems. 

Engineered quenches can provide a systematic way to study the relatively unexplored 
problem of inhomogeneous phase transitions. Using this technique we have shown how it 
is possible to invert the process and realise a homogeneous phase transition over a large 
spatial region of an harmonically trapped BEC. Although we were unable to drive the 
entire condensate through the critical point at the same time due to the breakdown of 
our approximation, the approach clearly establishes a path to simplify phase transitions in 
inhomogeneous systems. 

Finally we have addressed two specific issues of the KZ mechanism: the relation between 
causality and defect formation, and the critical moment at which the defect density is decided. 
For the former we have derived a threshold value for the velocity of a front at which defect 
formation is suppressed in our system. For the latter, using a temporally asymmetric quench, 
we have supplied additional evidence in support of the case made in Ref. [62], i.e. that 
the density of defects is decided after the transition for this particular realisation of the KZ 
mechanism. 
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Appendix A. Critical Coupling and Critical Exponents 



In this appendix we derive an expression for the critical coupling £l cl and critical exponent 
v for a uniform binary BEC in a ring trap with V(x) = (the definition of a uniform, 
homogeneous system). From Eq. ([T| follows that the time-independent Gross-Pitaevskii 
equation for ip^ is 

.2 Q2 

2 + 9n\^i\ 2 +S12|V'3-*| 2 Ipi ~ ntllpz-i. (A.l) 



2m dx 2 



Here we are interested in the behaviour of small perturbations around the mean-field solutions 
ipf of (A.l i. We decompose the wavefunction ipi into its mean field value ip® and the 
fluctuations (pi : ipi = -0° + 4>i- Without loss of generality we take the condensate mean-field 
functions ip® to be real. In order to simplify the equations further we work in the conditions 
where g = gu — 922 ^ .912 and Ni = N2 = N/2. The equality of the intraspecies 
interactions is a good approximation for the reality of mixtures of hyperfine states of 87 Rb, 
where the difference between the interaction constant is less than 5%. The equal number of 
particles can easily be realised experimentally by starting with a single component, before 



applying a 7r/2 pulse of the coupling. If we substitute the decomposition of ipi into Eq. (A.l 
and retain only terms linear in <j> we obtain the equation 



dx 2 



2 



(A.2) 
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where £b = ^/V^^rngp, 

S= ( ~ ■ "',',<, /,<>"' I (A.3) 




and p = N/L is the linear density. The diagonalisation of the matrix S results in a pair of 
eigenvalues 1^ = 2{ 1 + l/-\/A, 1 — 1/ a/A + Ml / gp} and a transformation matrix U suc h that 
A = USU^ 1 where A is a diagonal matrix. In the basis ((/)[, tfi' 2 ) T = U(<px, <p2) T , Eq. ( A.3 1 
takes the form 

We can then extract information about the miscible-immiscible phase transition by examining 
the effective correlation lengths £j = •v/^g/Ti- As expected, an uncoupled mixture (£1 = 0) 
has a critical point for A = 1. We are interested here in the value of the critical coupling 
which turns an immiscible condensate into a miscible one. From £2 it follows that when 
A < 1 the critical point is 

hn cI = gp(^-l). (A.5) 
v A 



It is also possible to quantify the critical exponent v by noticing that the effective correlation 
length scales as £2 = = / 1 e 1 1 ^ 2 which combined with Eq. ^ implies v — 1/2. 

Similar results for the critical coupling and critical exponent were obtained in Ref. Il49ll 
from the energy spectrum of the system. We quote here only the energy gap E gap oc 
hy/Q(t)[tt(i) — f2 cr ] of the spectrum in the mean-field approximation. The presence of a 
gapped spectrum indicates that the phase transition is continuous. The relaxation time is 
related to the energy gap through the relation r = h/E gap oc 1 / 1 e | 1 / 2 which implies the 
dynamical critical exponent is z = 1. 
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